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criterion  for  optimizing  the  parameters  of  the  2d  variational  (2dVar)  interpolation  of  high  frequency 
radar  (HFR)  data,  and  assessing  the  accuracy  of  the  surface  currents'  simulations  by  regional  models.  It  is 
shown  that  penalizing  the  magnitude  and  enforcing  smoothness  of  the  divergence  field  significantly 
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I.  Introduction 

In  recent  years  considerable  efforts  have  been  made  to  study 
predictability  of  the  trajectories  of  the  floating  material  at  the  sea 
surface  (Ozgbkmen  et  al.,  2001 ;  van  Sebille  et  aL.  2009;  Lumpkin 
and  Elipot.2010;  Huntley  et  al.,  2011).  This  problem  is  important  in 
many  applications,  including  search  and  rescue  missions,  accident 
localization  via  backward  tracking  of  the  debris,  or  optimization  of 
the  protective  missions  in  response  to  environmental  disasters. 

Since  numerical  modeling  became  the  major  forecast  tool  in 
(Keanography  a  laige  literature  appeared  on  the  Lagrangian 
predictability  (LP)  of  the  velocity  fields  produced  by  the  mcxlels. 
The  major  parameter,  characterizing  the  LP  of  a  model  is  the  root 
mean  square  (rms)  separation  e,  between  the  observed  and 
model  simulated  drifters  for  a  given  forecast  time  r.  A  few  years 
ago  Barron  et  al.  (2007)  conducted  a  comprehensive  LP  study  of 
the  Navy  Layered  Ocean  Model  by  comparing  drifter  tracks  in  30 
regions  of  the  World  Ocean  with  the  trajectories  simulated  by  the 
model.  The  estimated  values  of  e,  varied  in  the  range  between 
ei=10  25  km  for  T=lday  and  increased  to  50  150  km  for 
T  =  7  days.  Later,  Huntley  et  aL  (2011 )  studied  the  LP  of  the  Navy 
Coastal  Ocean  Model  (NCOM)  at  1/16°  resolution  in  the  South 
China  Sea  documenting  the  values  of  ei  =  17  km  and  e7=145  km. 
They  also  observed  that  coarsening  of  spatial  resolution  up  to  1/2° 
had  a  relatively  small  impact  on  the  LP.  In  contrast,  reduction  of 
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the  tempioral  resolution  to  12  20  h  significantly  reduced  the  LP 
due  to  poor  resolution  of  tidal  currents.  In  the  most  recent  study. 
Scott  et  al.  (2012)  have  shown  that  LP  of  the  model  output  could 
be  improved  by  taking  the  ensemble  average  of  the  drifter 
trajectories  produced  by  several  models.  In  particular,  the  viilue 
of  Cl  was  reduced  from  27  km  to  21  km  when  averaging  over  five 
models  of  the  equatorial  Atlantic  was  performed. 

In  general,  it  appears  that  a  typical  model  value  of  e-i  varies 
within  15  25  km  and  rarely  falls  below  10  km  (in  the  latter  case,  as 
a  rule,  the  background  flow  is  strong  and  exhibits  small  shear).  It  is 
also  notewortfty  that  LP  is  weakly  affected  by  data  assimilation 
(Scott  et  aU  2012),  unless  independent  velocity  observations  are 
taken  along  the  drifter  trajectories. 

In  this  respect,  development  of  the  coastal  HFR  networks 
(Harlan  et  al.,  2010)  can  be  considered  as  a  major  observational 
support  for  increasing  the  accuracy  of  the  Lagrangian  forecasts  in 
the  near  shore  regions.  HFR  observations  are  capable  of  delivering 
information  on  surface  currents  that  is  consistent  in  both  space  (1 
5  km)  and  time  (0.2  1  h)  resolutions  of  the  regional  circulation 
models.  Therefore,  in  recent  years  a  lot  of  studies  have  been 
devoted  to  LP  comparisons  of  the  surface  currents  derived  from 
HFR  observations  by  various  techniques.  As  a  few  examples. 
Ullman  et  al.  (2006)  investigated  the  skill  of  the  HFR  derived 
currents  in  predicting  drifter  trajectories  on  the  Atlantic  coast  of 
the  US  and  obtained  the  value  of  ei=5  7  km,  using  Lagrangian 
tracking  with  a  random  flight  model  to  simulate  the  effect  sub 
grid  turbulence  on  the  trajectories:  Kohut  et  al.  (2012)  employed 
8  drifters  to  assess  the  parameters  of  the  optimal  interpolation 
scheme  which  was  used  to  obtain  gridded  velocities  from  the  HFR 
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observations  off  the  New  Jersey  coast,  obtaining  the  value  of 
ei=8  km;  Shadden  et  aL  (2009)  have  shown  a  substantial  coher 
ence  between  the  drifter  trajectories  and  the  HFR  derived  Lyapu 
nov  exponents  in  the  Monterrey  Bay.  In  their  recent  work,  Barrick 
et  al.  (2012)  demonstrate  a  good  forecast  skill  of  a  simple  statistical 
model  of  the  surface  currents  derived  from  HFR  data  by  the  OMA 
technique  (Kaplan  and  Lekien,  2007). 

These  results  indicate  that  HFR  observations  have  a  better 
forecasting  capability  on  their  own  right,  and  may  substantially 
increase  the  LP  of  the  surface  velocity  fields  after  being  assimilated 
into  the  numerical  models. 

In  the  present  study  we  assess  LP  of  the  HFR  derived  currents 
during  the  Deep  Water  Horizon  (DWH)  oil  spill  (May  December 
2010).  In  particular,  we  explore  the  impacts  on  the  LP  of  the 
preliminary  gap  filling  and  of  suppressing  divergence  of  the  HFR 
currents  obtained  by  the  2d  variational  algorithm  (Yaremchuk  and 
Sentchev,  2009, 2011).  The  LP  is  used  as  a  benchmark  in  estimating 
the  overall  accuracy  of  HFR  observations  in  the  DWH  region  and 
optimization/tuning  the  2dVar  interpolation  scheme.  The  potential 
benefit  of  assimilating  HFR  data  into  regional  models  is  assessed 
by  comparing  the  LPs  of  the  HFR  derived  velocity  fields  with  the 
data  assimilative  NCOM  output 

The  paper  is  organized  as  follows.  In  the  next  section,  we 
delineate  the  study  area,  and  describe  the  data.  In  Section  3.  the 
methodologies  used  to  estimate  surface  currents  and  the  respec 
tive  Lagrangian  tr^ectories  are  outlined.  Then  in  Section  4  we 
compare  the  trajectory  forecast  skill  statistics  for  various  HFR 
interpolation  products  and  for  the  assimilative  NCOM  model  runs. 
We  conclude  with  a  summary  of  the  key  results  and  their 
implications  for  further  studies  of  the  DWH  region. 


2.  Observations  in  May-December  2010 
2.1.  HFR  observations 

Radial  sea  surface  currents  were  measured  by  three  CODAR 
SeaSonde  systems  operating  along  the  northern  Gulf  of  Mexico 
(NCOM)  coast  in  the  region  east  of  the  Mississippi  mouth.  The 
three  site  long  range  (-4.5  MHz)  system  provided  hourly  surface 
currents  at  6  km  radial  and  5°  azimuthal  resolution  over  the  region 
shown  in  Fig.  1.  Radial  velocities  observed  at  each  of  the  3  sites 
were  obtained  with  the  MUSIC  algorithm  (Schmidt,  1986)  using 


Fig.  1.  HFR  coverage  of  the  DWH  region.  Dots  show  locations  of  the  radial  velocity 
observations  from  the  three  radars  denoted  by  gray  triangles  Intensity  of  the  gray 
color  of  a  dot  is  proportional  to  the  total  observation  time  during  the  period 
between  7  May  and  31  December.  2010.  Shaded  areas  delineate  subdomains 
covered  by  HFR  observations  more  than  80%  of  time  (dark  gray),  65 — 80%  of  time 
(gray)  and  50 — 65%  of  time  (light  gray).  The  DWH  location  is  shown  by  the 
black  star. 


measured  antenna  patterns.  These  data  were  obtained  from  NOAA 
National  Data  Buoy  Center  HFR  archive. 

The  acquired  radial  velocity  data  were  characterized  by  inho 
mogeneous  coverage  both  in  space  and  time  (Fig.  2).  The  temporal 
patchiness  was  caused  by  multiple  reasons,  mostly  by  malfunc 
tions  of  the  radars  (e.g.,  in  the  beginning  of  August)  and  required 
filling  the  gaps  in  the  space  time  coverage  of  the  domain  by 
observations  prior  to  retrieving  the  gridded  velocity  fields  from 
the  data.  The  utilized  gap  filling  method  (Yaremchuk  and 
Sentchev.  2011)  is  similar  to  the  technique  used  in  preprocessing 
satellite  SST  data  that  are  often  obscured  by  clouds.  After  filling  of 
the  gaps,  hourly  velocity  vector  fields  were  reconstructed  using 
the  variational  interpolation  method  described  in  the  next  section. 

The  effective  depth  of  the  measured  surface  currents  depends 
on  the  HFR  frequency.  Under  the  assumption  of  a  linear  surface 
current  vertical  profile  the  effective  measurement  depth  is  pro 
portional  to  the  radar  wavelength  with  the  proportionality  coeffi 
cient  of  1  f&x  (Stewart  and  Joy.  1974).  Therefore,  the  effective  depth 
of  the  current  measurements  is  estimated  to  be  -2.6  m. 

Assessment  of  the  uncertainties  provided  by  the  CODAR  sys 
terns  was  made  under  the  assumption  that  the  respective  errors 
a  are  <5  correlated  in  space.  The  error  estimates  were  made  using 
statistical  analysis  of  the  radial  velocities  in  the  course  of  the  gap 
filling.  The  resulting  values  of  <r  were  computed  as  square  roots  of 
the  azimuthally  averaged  diagonal  elements  of  the  noise  covar 
iance  matrix  and  varied  from  4  cm/s  to  16  cm/s  for  the  ranges 
between  15  and  240  km.  These  numbers  are  consistent  with  the 
typical  error  estimates  of  the  long  range  CODAR  systems.  The 
values  of  a  were  used  in  the  definition  of  the  cost  function  weights 
in  the  2dVar  algorithm  for  retrieving  the  vector  fields  from  the 
observed  radial  components. 


2.2.  Drifter  trajectories 

We  used  the  observed  trajectories  of  satellite  tracked  drifters 
available  at  the  Naval  Oceanographic  Office  database  for  the 
8  months  (May  December)  of  2010.  The  database  contains  drifter 
data  from  many  sources  including  the  AOML  Drifter  Data  Assem 
bly  Center  (www.aomLnoaa.gov/phod/dac/dacdata.html  and  the 
US  Coast  Guard.  In  addition  we  used  quality  controlled  drifter 
trajectories  from  the  Ocean  Circulation  Group  at  the  University  of 
South  Florida.  These  data  included  drifters  that  were  deployed  in 
the  eastern  Gulf  of  Mexico  region  in  summer  2010  as  a  rapid 
response  to  the  Deepwater  Horizon  oil  spill  (Liu  et  al.  2011; 
Weisberg,  2011). 

Most  of  the  drifters  were  drouged  to  follow  currents  near  the 
12m  depth,  which  is  consistent  with  the  effective  depth  range  of 
the  currents  registered  by  the  HFRs.  The  drifters  have  been  tracked 
with  several  satellites,  with  1  h  time  between  the  fixes.  We  used 
the  quality  controlled  drifters  interpolated  to  1  hourly  positions  in 
the  time  interval  between  May  7,  and  December  31  of  2010.  In 
space,  the  domain  was  limited  to  the  area  of  HFR  coverage  (Fig.  1). 
All  the  drifter  tracks  considered  were  located  within  the  regular 
6  km  grid  (Fig.  3)  used  for  HFR  data  interpolation  (see  Section  3.1). 
Drifter  trajectories  observed  between  May  7  and  December  31. 
2010  in  the  study  area  are  plotted  in  Fig.  3. 

On  the  total,  37  drifters  (comprising  297  drifter  days)  were 
used  in  the  analysis.  It  should  be  noted  that  the  major  part  of  the 
trajectories  are  located  in  the  regions  with  relatively  poor  HFR 
coverage  (cf.  Figs.  1  and  3):  many  drifters  were  released  immedi 
ately  after  the  spill  in  the  vicinity  of  the  DWH  site  (Fig.  3),  which 
was  covered  by  coastal  radars  only  55%  of  the  time  (Fig.  1). 
Temporal  coverage  was  also  uneven:  75%  of  the  data  (226 
drifter  days)  were  acquired  in  the  period  between  May  7  and  July 
31.  2010. 
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Fig.  X  Drifter  tracks  for  the  S-month  period  of  the  study  (M^  7, 2010  to  Dec  31. 
2010).  The  6  km  interpolation  grid  for  the  HFR  data  is  shown  in  gr^.  The  DWH 
location  is  shown  by  the  black  star. 


23.  NCOM  model  output 

In  the  present  study  we  used  two  data  constrained  solutions  of 
the  NCOM  model  (Martia  2000;  Barron  et  al..  2006:  Martin  et  aL, 
2008).  The  employed  version  of  NCOM  uses  a  combined  a  z  vertical 
grid  with  32  levels  and  the  Arakawa  C  grid  in  the  horizontal  at  3  km 
resolution.  The  data  assimilation  scheme  (Cummings,  2005)  cur 
rently  uses  a  3d  variational  method,  which  allows  us  to  assimilate 
satellite  observations  of  sea  surface  height  and  temperature  and 
in  situ  temperature  and  salinity  data  from  various  sources. 

The  first  set  of  data  used  in  our  experiments  is  the  combined 
set  of  3dVar  NCOM  analyses  (computed  daily  at  0  UTC)  and  6  h 
forecasts  (at  6, 12  and  18  UTC)  generated  by  integrating  the  model 
from  the  analysis  times.  The  hourly  velocity  components  were 
obtained  by  linearly  interpolating  these  four  model  outputs  within 
the  24  h  window.  Since  the  6  km  HFR  interpolation  grid  (Fig.  3) 
was  configured  to  coincide  with  every  other  node  of  the  NCOM 
grid,  the  output  currents  used  for  Lagrangian  tracking  were 
obtained  by  [>icking  every  other  grid  point  value  from  the  original 
model  grid.  These  data  cover  the  whole  period  from  May  7  to 
December  31.  2010.  They  will  be  further  abbreviated  by  rNCOM 
(regional  NCOM). 

The  second  set  of  velocity  fields  is  the  ensemble  mean  of  the 
NCOM  runs  constrained  by  the  above  described  data  using  the 
ensemble  transfer  (ET)  technique  (Bishop  and  Toth.  1999).  The 
ensemble  contained  32  members,  generated  by  perturbing  atmo 
spheric  forcing  derived  from  the  Navy  atmospheric  model  using 
the  time  deformation  technique  (Wei  et  al..  2013).  Additionally, 
both  vertical  and  horizontal  mixing  parameters  were  perturbed  to 
represent  model  uncertainties  from  NCOM  sub  grid  parameteriza 
tions.  The  mixing  parameter  perturbations  were  Gaussian  with  the 
means  of  0.125, 17.5  and  variances  of  0.01875  and  0.625  for  the 
Smagorinsky  and  Mellor  Yamada  coefficients  respectively.  More 


details  on  the  ET  implementation  with  NCOM  can  be  found  in  Wei 
et  aL  (2013).  The  set  of  the  ensemble  analyses/forecasts  (denoted 
eNCOM  further  below)  was  available  on  the  same  space  time  grid 
as  rNCOM,  but  for  the  shorter  period  of  time  from  May  7  to  July  25, 
2010.  Similar  technique  was  used  to  project  the  eNCOM  near 
surface  velocities  on  the  HFR  interpolation  grid.  Hourly  velocities 
in  both  model  outputs  were  taken  from  three  different  vertical 
levels:  0  m,  25  m  and  7.5  m. 


3.  Methodology 
3.1.  Processing  HFR  data 


Since  HFRs  measure  projections  of  the  surface  velocity  vectors 
on  the  directions  r  of  the  radar  beams,  an  algorithm  is  needed  to 
reconstruct  the  velocity  field  from  such  observations.  In  the 
present  study  we  use  a  combination  of  the  statistical  gap  filling 
technique  with  the  2dVar  interpolation  algorithm  (Yaremchuk  and 
Sentchev,  2011)  for  this  purpose.  The  approach  is  esp>ecially 
suitable  for  gappy  data  (Fig.  2)  with  inhomogeneously  distributed 
observation  points  and  allows  to  directly  control  divergence  and 
vorticity  of  the  interpolated  field  u. 

The  optimal  estimate  of  u  is  determined  by  maximizing  its 
likelihood  expressed  in  terms  of  the  Gaussian  probability  density 
function  7Yu)~exp[-y(u)].  The  argument  of  the  exponent  (the  cost 
function)  is  quadratic  in  u  and  consists  of  two  terms:  J=Jd+  Jr- 
The  first  term  Ja  measures  the  distance  between  the  radial 
velocity  uj  observed  at  the  kth  point  and  the  corresponding 
components  of  the  interpolated  field  u  at  the  observation  points: 


i  crj-2[(U.rk)  1^2 

k  I 


(1) 


Because  of  Gaussianity  the  distance  is  quadratic  in  u  and  scaled  at 
the  observation  points  by  the  corresponding  measurement  error 
variances  a^.  The  regularization  term  Jr  is  introduced  to  penalize 
higher  spatial  derivatives  of  u: 

Jr=-^  f  [Wu(Au)^  +  Wrf(A  div  uf  +  WciA  curl  u)^]  dx  dy  (2) 

where  divu  =  5xU+SyV,  curl  u  =  o!»v-ayU  are  the  divergence  and 
vorticity  respectively,  A=dxx  +  dyyis  the  Lapladan  operator  and  A 
is  the  area  of  interpolation  domain  shown  in  Fig.  3.  The  differential 
operators  were  approximated  by  central  finite  differences  on  the 
regular  interpolation  grid  with  a  step  <5x  =  6  km  equal  to  the  radial 
resolution  Ir  of  the  HFR  data. 

The  weights  Wu,Wd  and  Wc  have  the  meaning  of  the  inverse 
error  variances  of  the  corresponding  fields  and  allow  to  control 
their  magnitude  in  the  interpolation  pattern.  In  the  experiments 
described  in  the  next  section,  the  cost  function  weights  were 
derived  from  HFR  data  statistics  and  then  fine  tuned  to  maximize 
the  predictability  of  the  drifter  tr^ectories. 
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Another  important  aspect  of  the  HFR  data  processing  is  gap 
filling.  The  back  scattered  HFR  signals  suffer  from  distortions  of 
artificial  and  natural  origin,  such  as  atmospheric  conditions,  sea 
surface  roughness,  sea  traffic  and  mallunctions  in  radar  operation. 
As  a  consequence,  estimates  of  the  radial  velocities  extracted  from 
the  Doppler  shifts  of  the  HFR  signals  become  unusable,  resulting  in 
numerous  gaps  in  spatial  coverage,  which  was  the  case  with  the 
data  we  used  (Fig.  2).  As  it  is  seen  from  Figs.  1  and  2.  most  of  the 
gaps  have  relatively  short  life  time  and  are  more  frequent  at  long 
ranges.  At  times,  a  large  amount  of  data  is  lost  due  to  improper 
operation  of  a  single  radar  when  the  radial  velocities  are  still 
observed  at  nearby  locations  by  the  two  remaining  HFRs. 

In  this  case  the  gap  filling  technique  based  on  the  EOF  decom 
position  of  the  sample  covariance  matrix  (e.g..  Beckers  and  Rixen. 
2003)  is  quite  efficient  In  the  present  study  we  employed  a  similar 
method  adapted  for  HFR  observations  by  Yaremchuk  and  Sentchev 
(2011 ).  The  method  usually  works  well  for  the  observation  loca 
tions  where  gaps  occupy  less  than  f~20  30%  of  the  total  observa 
tion  time.  In  our  case  this  condition  for  the  gap  concentration 
parameter  f  was  satisfied  for  approximately  60%  of  the  points 
where  radial  velocities  were  observed. 

The  gap  filling  procedure  was  performed  in  several  steps.  First  a 
set  of  cross  validation  (CV)  points  randomly  distributed  across  the 

8  month  observation  period  was  removed  from  the  data.  On  the 
average.  1%  of  the  data  were  set  aside  from  each  hourly  observation 
of  the  radial  velocities.  Second,  iterative  processes  of  estimating  the 
covariance  matrix  and  the  gap  filling  were  performed  for  different 
numbers  M  of  the  gap  filling  EOFs.  Third,  the  optimal  value  of  M  was 
selected  by  minimizing  the  gap  filling  error  eg  which  was  assessed 
using  the  temporarily  removed  observations  in  the  CV  points.  Finally, 
the  CV  points  were  added  to  the  data  set  and  the  remaining  gaps 
were  filled  using  the  optimal  number  of  EOFs 

It  is  noteworthy  that  the  optimal  value  of  M  allows  for  separa 
tion  of  signal  from  noise  in  the  data  space  (e.g..  Yaremchuk  and 
Sentchev.  2011):  spatial  structures,  described  by  M  EOFs  with 
largest  amplitudes  can  be  attributed  to  the  “resolved  signal" 
whereas  variability  described  by  the  rest  of  the  EOFs  can  be  treated 
as  noise.  In  our  experiments  the  optimal  value  of  M  varied  in  the 
range  between  25  and  64.  accounting  for  60  80%  variabQity  of  the 
radial  velocities.  The  noise  covariance  matrix  C  had  negligible 
correlations  between  the  radial  velocity  errors,  so  we  used  the 
inverse  diagonal  elements  of  C  for  weighting  the  misfits  between 
the  observed  and  interpolated  radial  velocities  (Eq.  (1)). 

In  the  course  of  numerical  experiments  the  gap  filling  was 
performed  in  three  sets  of  points  characterized  by  three  values  of 
4  =  50, 30  and  1 5%.  The  filled  points  contributed  respectively  by  28, 

9  and  1.3%  to  the  total  number  (6  324  530)  of  the  observed  radial 
velocities.  After  normalization  by  the  rms  variance  of  the  observed 
radial  velocities  the  gap  filling  errors  varied  between  056  and 
0.64  for  the  different  values  of  4. 


Fig.  4.  Examples  of  the  mean  virtual  drifter  tr^ecdaries  (gray)  used  for  assessing 
the  LP  skill  and  the  corresponding  real  drifter  tracks  (black).  The  clouds  of  light  gr^ 
dots  demonstrate  the  spread  of  the  virtual  tr«geaories  around  the  mean  at  the  end 
point  Launch  points  of  the  virtual  drifters  are  shown  by  black  circles.  Numbers 
show  the  integration  time  in  days. 


(e.g.,  Ullman  et  al..  2006)  with  the  effective  diffusion  coefficient 
of  1.5  X  10®cm^/s.  a  typical  value  of  the  Smagorinsl^r  diffusion 
coefficient  from  the  NODM  model.  This  approach  provides  a  con 
sistent  Lagrangian  tracking  techniques  for  the  velocity  fields  obtained 
from  both  HFR  observations  and  the  assimilative  models  (Fig.  4). 

The  prediction  of  the  particle  trajectory  was  done  by  integrating 
Eq.  (3)  with  a  fourth  order  Runge  Kutta  scheme.  At  hourly  times  r 
of  each  integration  of  (3).  the  prediction  error  e,  was  computed  as 
the  distance  between  the  location  Xc  of  the  centroid  of  the  particle 
cluster  and  the  actual  drifter  position  :  e,  =  |Xc(t)-Xj(t)|.  Predic 
tion  errors  were  then  averaged  0  over  all  the  available  segments, 
whose  number  n  varied  between  71  (t=  1  h)  and  20  (t  =  9  days)L 
Since  the  segment  length  (4.5  days)  was  chosen  as  the  mean 
decorrelation  time  scale  of  the  HFR  derived  (4.3  days)  and  model 
(4.7  days)  velocity  fields,  prediction  error  estimates  for  each  seg 
ment  were  assumed  to  be  statistically  independent.  Error  bars  for 
(e>  were  computed  as  the  70%  confidence  intervals  derived  from  the 

distributions  with  2n  degrees  of  freedom. 

In  addition  to  the  virtual  particle  trajectories,  the  mean 
separation  <|Xd(T)-Xd(0)|>  between  the  actual  drifter  position  at 
time  r  and  its  initial  position  Xd(0)  at  the  segment  was  also 
computed.  The  prediction  skill  y  was  defined  by 


<er) 

<|Xd(T)  x;,(0)|> 


(4) 


Thus,  the  skill  was  estimated  relative  to  the  “position  persis 
tence  scenario"  which  represents  a  search  based  on  the  last  known 
position  Xd(0)  in  the  complete  absence  of  information  on  the 
surface  currents. 


3.2.  Computing  trajectories  and  their  predictability 

Drifter  tracks  were  used  to  evaluate  errors  of  the  partide 
trajectories  corresponding  to  the  velodty  fields  u  that  were  obtained 
either  from  HFR  observations  or  from  the  model  Each  drifter 
trajectory  was  divided  into  45  day  segments  and  the  actual  drifter 
positions  at  the  beginning  of  the  segments  were  taken  as  the  initial 
condition  for  the  virtual  particle  release.  1000  particles  were  released 
within  each  segment  into  the  2d  velodty  field  described  by 

dx 

^=H(X.t)+ U'(X.t)  (3) 

where  X  is  the  position  of  the  partide  and  u'  is  the  stochastic 
contribution  of  the  sub  grid  scale  variability  iind  the  measurement 
uncertainty.  The  model  used  to  determine  u'  is  a  random  walk  model 


4.  Results 

Computations  reported  in  this  section  had  two  objectives.  The 
first  one  was  to  explore  how  preprocessing  of  the  HFR  data  affects 
the  LP  skill  of  the  reconstructed  velocity  fields.  The  second  one 
was  to  compare  in  terms  of  the  LP  skill  the  HFR  derived  velocities 
with  the  available  operational  data  assimilation  products,  uncon 
strained  by  the  surface  velocity  data. 

4.1.  Prediction  skill  of  the  HFR  derived  currents 

4.1.1.  Impact  of  the  prior  statistical  assumptions 

As  it  has  been  mentioned  in  Section  3.1,  result  of  2dVar 
interpolation  crucially  depends  on  the  statistics  of  the  interpolated 
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fields,  which  is  expressed  in  terms  of  the  inverse  error  variances 
a  Wu.Wj  and  Wc  of  the  cost  function. 

The  interpolation  error  variances  which  measure  the  misfit 
between  the  observed  radial  velocities  and  their  counterparts 
derived  from  the  interpolated  pattern  (Eq.  (1))  were  estimated  as 
the  diagonal  values  of  the  noise  covariance  matrix  (Section  3.1). 
Their  values  varied  between  16  cm^/s^  in  the  close  proximity  from 
the  radars  to  about  150  250cm^/s^  at  far  ranges  (150  240  km). 
The  regularization  weights  (Eq.  (2))  were  first  estimated  by 
considering  the  typical  spatial  scales  of  the  reconstructed  velocity 
field  and  then  fine  tuned  to  obtain  the  best  skill  in  predicting  real 
drifter  trajectories. 

As  a  rough  estimate  of  Wu,  we  assumed  that  the  smallest 
spatial  scale  /  is  twice  the  range  discretization  l~2lr  =  12  km  of  the 
radial  velocities.  This  assumption  translates  into  the  estimate 
/~(<r^/Wu)'/‘'  =  2/r,  or  Wu~(T^/f/16.  Rne  tuning  of  this  weight 
against  the  drifter  data  (with  Wc=  Wj  =  0)  has  shown  that  the 
best  LP  skill  is  achieved  at  Wu  =  0.05<r^t^,  a  value,  corresponding  to 
the  cut  off  scale  2.1 5/r.  The  resulting  mean  interpolation  error 
(e)  =  (<yKJj)  was  0.32,  fairly  consistent  with  the  estimated  noise 
level  in  the  radial  velocity  data. 

After  fixing  the  optimal  value  of  <e>,  a  large  set  of  2dVar 
interpolation  experiments  has  been  conducted  to  optimize  the 
values  of  Wc  and  Wd.  These  experiments  have  shown  that  (y>  (the 
value  of  fr  averaged  over  the  forecast  times  between  0  and  9  days) 
is  weakly  sensitive  to  Wc  and  exhibits  a  slight  deaease  as  the  value 
of  Wc  approaches  zero.  For  this  reason  the  corresponding  term  was 
eliminated  from  the  cost  function. 

In  contrast  to  Wc,  varying  Wd  which  suppresses  the  roughness 
and  magnitude  of  the  divergence,  had  a  significant  impact  on  the 
LP  skill  at  the  forecast  times  exceeding  3  4  days  (Fig.  5).  The 
optimal  value  of  Wd  was  found  to  be  8  x  10  ^  cm''  s^,  and  resulted 
in  the  2  3  smaller  rms  magnitude  of  the  divergence  compared  to 
the  curl  of  the  surface  currents.  The  better  LP  skill  at  larger  forecast 
times  can  be  partly  explained  by  the  fact  that  suppressing  the 
divergence  tends  to  increase  the  accuracy  in  reconstructing  the 
geostrophic  component  of  the  flow,  which  is  characterized  by 
slower  variability  than  tidal  and  wave  induced  components  of  the 
surface  currents. 

Rg.  5  provides  an  indication  that  despite  poorer  Lagrangian 
statistics  at  longer  time  scales,  the  effect  appears  to  be  statistically 
significant:  the  error  reduction  observed  at  t>3  days  persists 
throughout  the  entire  forecast  time  interval. 

In  terms  of  the  absolute  value,  the  LP  error  grows  from  12 
20  km  at  the  forecast  times  of  1  2  days  to  35  50  km  at  t=5  9 
days.  The  short  term  forecast  errors  are  consistent  with  the  recent 
LP  study  of  Kohut  et  al.  (2012)  who  validated  HFR  derived  currents 
against  six  drifters  in  the  region  of  the  Mid  Atlantic  Bight  and 
obtained  the  values  of  8  14  km  for  the  forecast  times  of  1  2  days. 


Fig.  5.  Difference  between  the  forecast  errors  of  the  2dVar  interpolated  HFR 
fields  obtained  with  0  and  with  various  values  of  the  gap  concentratbn 
parameter  70%  confidence  limit  for  the  difference  is  shown  by  the  thick  dotted 
line.  Positive  values  correspond  to  smaller  forecast  error  of  the  divergence- 
penalizing  interpolation  (W^^). 


Fig.  6.  Drifter  forecast  skill  y,  as  a  function  of  the  forecast  time  r  for  the  results  of 
2dVar  interpolation  of  HFR  data  with  various  values  of  the  gap  concentration 
parameter  (  for  the  case 


When  referenced  to  the  persistence  errors  (Eq.  (4)),  these  numbers 
correspond  to  the  values  of  z,  around  0.65  0.7,  that  is  25  30% 
better  than  those  shown  in  Fig.  6.  Taking  into  the  account  more 
complex  (less  deterministic)  structure  of  the  flow,  poorer  quality 
of  the  HFR  data  (Fig.  2)  and  confinement  of  most  of  the  drifters  to 
the  outer  regions  of  the  HFR  footprint  (Fig.  1),  the  LP  skill  of  the 
2dVar  algorithm  can  be  considered  as  satisfactory. 


4.1.2.  Impact  of  the  gap  filling 

Rg.  6  shows  that  the  forecast  skill  gradually  degrades  with  the 
decrease  of  the  parameter  ^  which  defines  the  upprer  limit  of  the 
gap  percentage  for  filling  an  observation  point  with  a  value 
derived  from  statistical  analysis. 

Because  of  the  extremely  intermittent  time  coverage  (Rg.  2) 
only  1.5%  of  the  observation  points  were  characterized  by  less  than 
15%  of  gaps  (^<0.15)  and  all  of  these  points  were  located  at 
relatively  short  (20  50  km)  distances  from  the  radars,  whereas 
most  (75%)  of  the  drifter  positions  were  observed  at  far  ranges 
(120  180  km)  where  the  gap  concentration  was  fairly  high 
(f-0.4  0.7)  and  the  accuracy  of  HFR  observations  is  low. 

Table  1  shows  improvement  (%)  of  the  key  parameters  of  the 
velocity  field  normalized  by  their  values  obtained  without  pre 
liminary  gap  filling.  Apart  from  the  forecast  time  averaged  values 
of  Y,  and  e,  we  computed  the  trace  of  the  correlation  matrix  C  and 
the  relative  error  between  the  HFR  derived  velocities  Ur  and  the 
Eulerian  velocity  estimates  Vj  obtained  from  finite  differentiation 
of  the  drifter  positions  in  time: 


trC=C(Ur,Ud)+C(Vr,Vd); 


er=5 


l(Ur  Udf  +  (Vr  Vaf 


uj  +  v? 


Here  c  stands  for  the  correlation  coefficient  and  overline 
denotes  averaging  over  all  of  the  7117  Eulerian  velocities  estimated 
from  297  drifter  days  of  observations.  Zonal  Ua  and  meridional 
drifter  velocity  components  were  assessed  every  hour  using 
central  differences. 

As  it  is  seen  from  the  Table,  persistent  improvement  is 
observed  only  for  the  case  of  ^  =  0.15,  and  its  magnitude  (1.2 
1.6%)  approximately  corresponds  to  the  relative  inaease  in  the 
number  of  the  radial  velocities  subject  to  interpmlation  (1.3%)  after 
the  gap  filling.  With  the  increase  of  4,  the  number  of  additional 
(gap  filled)  observations  grows  dramatically,  but  they  do  not  result 
in  further  improvement  of  the  parameters.  Instead,  the  gap  filling 
“observations"  tend  to  decrease  the  proximity  between  the  HFR 
derived  velocity  fields  and  drifter  velocities. 

This  can  be  explained  by  the  above  mentioned  mismatch 
between  the  spatial  distribution  of  the  drifter  trqectories  and 
the  density  if  HFR  observations.  With  a  more  homogeneous 
temporal  HFR  and  spatial  drifter  coverages  than  shown  in 
Rgs.  2  and  3  one  may  expect  better  performance  of  the  gap 
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Table  1 

Changes  in  the  error  e^  correlation  ti(CX  and  the  mean  LP  skills  O'),  <e)  averaged 
over  the  forecast  times  after  the  gap-filling  performed  with  ^rious  threshold  gap 
concentrations  Numbers  (%)  are  normalized  by  the  results  obtained  for  the 
veloci^  fields  retrieved  by  the  2dVar  without  preliminary  gap-filling.  Positive 
values  denote  improvement  (decrease  of  er.<Cr>»0',>  and  increase  of  ti(C). 


Wd 

i 

Cr 

ti(C) 

O'.) 

<o 

0 

0.15 

13 

0.6 

0.4 

13 

1.9 

030 

9.1 

03 

33 

-43 

-2.8 

0.50 

2a7 

-2.4 

-1.1 

-7.8 

-7.6 

8 

0.15 

13 

1.1 

0.8 

13 

2.0 

0.30 

9.1 

03 

1.6 

-3.7 

-1.9 

0.50 

2S.7 

-2.5 

3.7 

-5.7 

-4.8 

filling  technique,  possibly  up  to  the  values  of  f~0.3.  This  requires 
additional  experimentation  with  the  appropriate  data  sets. 


42.  Comparison  with  the  NCOM  output 

Parameters  listed  in  Table  1  were  also  computed  for  the  available 
surface  and  subsurface  (0.  2  and  75  m)  velocity  fields  of  the  two 
NCOM  data  assimilation  runs  described  in  Section  23.  The  2  m 
velocities  demonstrated  the  best  values  of  the  parameters.  However, 
the  HFR  derived  velocity  fields  were  still  25  40%  better  both  in  terms 
of  the  LP  (Rg.  7)  and  discrepancies  between  the  Eulerian  velocities. 

The  NCOM  LP  error  curves  in  Fig.  7  compare  quite  favorably 
with  the  results  of  Price  et  al.  (2006)  who  validated  a  combination 
of  surface  velocities  from  the  ORSA  model  (Smith  et  al..  1982)  and 
ECMWF  wind  induced  currents  in  the  period  of  1997  1999  in 
approximately  the  same  region.  The  diagnosed  average  discrepan 
cies  were  found  to  be  78  km  (r  =  3  days).  229  km  (t  =  10  days)  and 
400  500  km  (r  =  20  30  days).  It  should  be  noted,  however,  that 
the  domain  considered  by  Price  et  al.  (2006)  extended  farther 
offshore  and  included  regions  dominated  by  mesoscale  eddies  and 
currents  of  the  open  Gulf. 

Compared  to  the  HFR  derived  velocities,  velocity  error  Cr  for  the 
NCOM  fields  was  approximately  10  15%  larger,  and  the  correlation 
coefficient  was  smaller  by  30%. 

In  general,  both  data  driven  solutions  in  Fig.  7  demonstrate 
similar  skills,  except  for  the  forecast  times  exceeding  7  days,  where 
the  ensemble  solution  shows  15  20  km  larger  disaepancy  with 
the  drifters.  This  could  be  partly  attributed  to  somewhat  poorer 
statistics  of  eNCOM,  which  was  avaQable  only  in  the  May  June 
2010  compared  to  the  full  8  month  averaging  period  of  rNCOM. 

A  substantially  better  LP  skill  of  the  HFR  derived  currents  is  not 
surprising  because  the  considered  models  are  not  constrained  by 
surface  velocity  data  and  also  do  not  account  for  processes  affecting 
the  transFX)rt  by  near  surface  currents  such  as  Stokes  drift  and  small 
scale  turbulent  wind  wave  interactions.  Similarly,  the  NCOM  model 
has  a  limited  skill  in  accurate  simulation  of  the  spatial  stmcture  of  tidal 
and  inertial  oscillations  which  have  profound  signatures  in  the 
considered  regioa  ^nd  sometimes  dominate  geostrophic  and  wind 
driven  components  of  the  surface  flow.  In  this  sense  the  considered 
coastal  region  appears  to  be  more  challenging  from  the  LP  viewpoint, 
as  it  is  often  dominated  by  smaller  scale  processes  neither  well 
resolved  by  the  model  nor  manifested  in  the  wind  forcing.  As  a  result, 
model  generated  tr^ectories  have  littie  or  no  correlation  with  the 
drifter  trajectories.  In  fact,  the  typical  values  of  the  LP  skill  y  for  both 
NCOM  runs  were  dose  to  unity,  sometimes  reaching  13  15  (cf.  Fig  6) 
meaning  that  the  models,  on  the  average,  do  not  have  any  LP  skill 
compared  to  the  null  assumption  of  zero  surface  velodty.  This  result 
underlines  the  potential  benefits  of  assimilating  HFR  observations  into 
the  numerical  models  of  coastal  regions. 


Fig.  7.  Drifter  forecast  error  e,  for  the  results  of  2dVar  interpolation  of  HFR  data 
with  enforced  smoothness  in  the  divergence  field  (black  curve)  and  for  the  NCOM 
model  solutions  shown  by  thick  lines 


5.  Discussion  and  conclusions 

A  set  of  37  drifter  trajectories  was  used  to  evaluate  the 
Lagrangian  predictability  in  the  Deep  Water  Horizon  (DWH)  oil 
spill  region  based  on  local  HF  radar  data  and  two  data  assimilating 
regional  model  runs.  The  LP  was  used  as  a  criterion  to  search  for 
the  best  parameters  of  the  2dVar  algorithm  retrieving  the  surface 
velodty  fields  from  the  radial  velodties  recorded  by  the  HF  radars. 

It  is  shown  that  penalizing  the  magnitude  and  enfordng 
smoothness  of  the  divergence  field  tends  to  increase  LP  by  2  6% 
at  the  forecast  times  of  3  9  days  while  preserving  it  at  the  smaller 
forecast  times.  Applying  preliminary  gap  filling  technique  based 
on  the  analysis  of  spatial  correlations  of  the  radial  velocities  adds 
an  extra  1  2%  to  the  LP  and  improves  correlations  between  the 
HFR  derived  velocities  and  the  respective  Eulerian  velocities 
derived  from  drifters. 

The  analyzed  drifter  set  was  not  conforming  with  the  HFR 
coverage  pattern  as  most  of  the  drifters  were  released  in  the  DWH 
region  which  is  located  at  the  far  ranges  (150  200  km)  off  the 
radar  sites.  This  partially  explains  a  relatively  moderate  improve 
ment  in  the  LP  caused  by  gap  filling  of  the  observation  points 
covered  more  than  85%  of  the  time  and  failure  to  obtain  an 
improvement  if  much  more  numerous  observation  points  with 
gap  concentrations  less  than  30  50%  were  taken  into  account.  We 
also  assume  that  the  divergence  related  improvement  of  the  LP  at 
the  forecast  times  of  3  9  days  could  be  more  visible  with  a  more 
homogeneous  coverage  of  the  domain  by  drifter  trajectories. 

Comparing  the  LP  skills  of  the  HFR  derived  velodties  with 
surface  currents  taken  from  the  data  constrained  NCOM  solutions 
have  shown  25  40%  better  performance  of  the  observational 
product  The  result  is  not  surprising  if  previous  model/drifter 
diagnostic  studies  (e.g..  Price  et  aU  2006:  Barron  et  al..  2007)  are 
compared  with  validations  of  the  HFR  derived  velodties  against 
drifters  (e.g.,  Ullman  et  al.,  2006;  Huntley  et  al,  2011 ;  Kohut  et  al, 
2012).  We  assume  that  such  a  large  gap  in  perfomaance  is  due  to 
the  fact  that  most  of  the  tested  models  use  relatively  simple 
description  of  the  processes  responsible  for  transport  in  the  top 
2  m  of  the  water  column,  which  is  strongly  affected  by  the  waves/ 
swell  and  other  upper  layer  processes  controlled  by  the  momen 
turn  transfer  from  the  atmosphere.  In  that  respect,  assimilation  of 
the  HFR  data  into  numerical  models  is  the  most  straightforward 
way  to  improve  the  situation.  Unfortunately,  to  the  best  of  our 
knowledge,  there  has  been  no  assessment  of  the  HFR  data  impact 
on  the  LP  skill  in  the  known  assimilation  studies  involving  radar 
observations  (e.g,  Paduan  and  Shulman,  2004;  Hoteit  et  al.,  2009: 
Barth  et  aL,  2010;  Yu  et  al.,  2012). 

Comparison  of  the  model  runs  driven  by  different  data  assim 
ilation  algorithms  has  shown  virtually  no  difference  between  the 
RT  and  3dVar  methods.  We  attribute  the  departure  of  the  eNCOM 
curve  from  rNCXDM  at  t  >  7  d  in  Rg.  7  to  relatively  poor  statistics  of 
the  2.5  month  eNCOM  period  characterized  by  an  average  duration 
of  a  drifter  trajectory  of  16  days.  Given  this  limitation,  the  result 
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still  appears  to  be  inconsistent  with  the  conclusions  of  Scott  et  al. 
(2012)  who  have  shown  that  using  a  multi  model  ensemble  mean 
provides  a  significant  (20%)  improvement  of  the  LP  in  the  Equa 
torial  Atlantic.  The  inconsistency  could  be  caused  by  several 
reasons,  such  as  better  statistics  (-3000  vs  297  drifter  days), 
negligible  tides  and  much  coarser  (20  30  km  vs  3  km)  resolution 
in  the  LP  computations  of  Scott  et  al.  (2012).  It  should  be  also 
noted  that  the  LP  skill  demonstrated  by  NQ3M  in  the  DWH  region 
(ei=20km,  07=67  km.  Fig.  7)  is  still  better  than  the  best 
ensemble  mean  estimates  (e,=21  km,  67=122  km)  of  Scott  et  al. 
(2012)  in  the  Equatorial  Atlantic. 

Our  results  also  show  that  substantial  LP  improvements  can  be 
achieved  by  f>reprocessing  HER  data  with  a  combination  of  the 
gap  filling  technique  and  2dVar  interpolation  tuned  to  suppress 
the  divergence  field.  Apart  from  a  much  (20  30%)  better  LP 
compared  to  the  considered  data  driven  model  runs,  interpolated 
HER  radial  velocities  are  less  noisy  and  continuous  in  space  and 
time.  These  properties  make  them  attractive  for  assimilation  into 
numerical  models  as  the  interpolated  EIER  velocities  are  less  prone 
to  spurious  modes  generated  by  the  raw  radial  velocities  often 
characterized  by  highly  intermittent  patterns,  especially  at  the 
outer  regions  of  the  HER  footprint.  The  issue  of  efficiency  of 
assimilating  either  raw  radial  velocities  or  their  projections  on 
the  model  grid  remains  an  open  question,  which  requires  further 
research. 
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